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A new approach to study the nonhnear phase of gravitational clustering in an expanding uni- 
verse is explored. This approach is based on an integro-differential equation for the evolution of 
the gravitational potential in the Fourier space which is obtained by using a physically motivated 
closure condition. I show how this equation allows one to understand several aspects of nonlinear 
gravitational clustering and provides insight in to the transfer of power from one scale to another 
through nonlinear mode coupling. 



I. INTRODUCTION 

There is considerable amount of observational evidence to suggest that about 25 per cent of the energy density in 
the universe is contributed by self gravitating system of dark matter particles. The smooth, average, energy density 
of these particles contributes to the expansion of the universe while any small deviation from the homogeneous energy 
density will lead to gravitational clustering. One of the central problems in cosmology is to describe the non linear 
phases of this gravitational clustering starting from a initial spectrum of density fluctuations. It is often enough 
(and necessary) to use a statistical description and relate different statistical indicators (like the power spectra, nth 
order correlation functions etc.) of the resulting density distribution to the statistical parameters (usually the power 
spectrum) of the initial distribution. The relevant scales at which gravitational clustering is non linear are less than 
about 10 Mpc (where 1 Mpc « 3 x 10^^ cm is the typical separation between galaxies in the universe) while the 
expansion of the universe has a characteristic scale of about few thousand Mpc. Hence, non linear gravitational 
clustering in an expanding universe can be adequately described by Newtonian gravity provided the rescaling of 
lengths due to the background expansion is taken into account. 

As to be expected, cosmological expansion introduces several new factors into the problem as compared to the study 
of statistical mechanics of isolated gravitating systems. (For a general review of statistical mechanics of gravitating 
systems, see 0. For a sample of different approaches, see and the references cited therein. Review of gravitational 
clustering in expanding background is available in several textbooks in cosmology 0,3.) (a) The problem has now 
become time dependent and it will be pointless to look for equilibrium solutions in the conventional sense of the word, 
(b) On the other hand, the expansion of the universe has a civilizing influence on the particles and acts counter to 
the tendency of gravity to make systems unstable, (c) In any small local region of the universe, one would assume 
that the conclusions describing a finite gravitating system will still hold true approximately. In that case, particles in 
any small sub region will be driven towards configurations of local extrema of entropy (say, isothermal spheres) and 
towards global maxima of entropy (say, core- halo configurations). It is not clear how these effects and cosmological 
expansion interact with each other at intermediate length scales. 

Though this problem can be tackled in a 'practical' manner using high resolution numerical simulations, such an 
approach hides the physical principles which govern the behaviour of the system. To understand the physics, it 
is necessary to attack the problem from several directions using analytic and semi analytic methods. Several such 
attempts exist in the literature based on Zeldovich(like) approximations 5], path integral and perturbative techniques 
0, nonlinear scaling relations "T] and many others. In spite of all these it is probably fair to say that we still do not 
have a clear analytic grasp of this problem, mainly because each of these approximations have different domains of 
validity and do not arise from a central paradigm. 

The purpose of this paper is to attack the problem from a different angle, which has not received much attention 
in the past. The approach begins from the dynamical equation for the the density contrast in the Fourier space and 
casts it as an integro-differential equation. This equation is known in the literature (see, e.g. Q) but has received 
very little attention because it is not 'closed' mathematically; that is, it involves variables which are not natural to 
the formalism and thus further progress is difficult. I will, however, argue that there exists a natural closure condition 
for this equation based on Zeldovich approximation thereby allowing us to write down a closed integro-differential 
equation for the gravitational potential in the Fourier space. It turns out that this equation can form the basis for 
several further investigations. In fact, one purpose of this paper — which is probably somewhat more pedagogical 
than is usual — is to draw the attention of the community to this approach and encourage further investigations 
based on this equation and ansatz. 

Given the validity of the ansatz, several conclusions follow in a straight forward manner. For example, I will describe 
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how different aspects of non linear gravitational clustering (which are well known in the literature) can all be obtained 
by taking suitable limits of this equation. What is more important, the structure of this integro-differential equation 
allows one to understand the key issue of transfer of power from one scale to another in non linear gravitational 
clustering. If the initial power spectrum is sharply peaked in Fourier space, in a shell of radius |k| — ko, then the non 
linear gravitational clustering allows for both cascading (to smaller length scales) and inverse cascading (to larger 
length scales) of the power during the evolution. The formalism developed here shows that the inverse cascading 
leads to the well known fc^ tail at small wave numbers (that is, at large spatial scales). Again, while this result is 
known, the derivation given here is new and appears to be simple, straightforward and does not require any other 
extra ad-hoc assumption. The cascading of power, on the other hand, proceeds broadly through the generation of 
harmonics at 2fco,4fco etc. with the mode couplin gq uickly leading to a universal power spectrum. These results 
agree with numerical simulations done in the past Finally, it turns out that one can obtain scale free solution 
to the integro-differential equation and analyse its asymptotic properties completely. This study shows some striking 
similarities between nonlinear gravitational clustering of coUisionless particles in an expanding universe and standard 
fluid turbulence. This application as well as its verification in numerical simulations will be described in a separate 
paper 



II. NONLINEAR GRAVITATIONAL CLUSTERING 



The expansion of the universe sets a natural length scale (called the Hubble radius) dn = c{a/a)~^ which is about 
4000 Mpc in the current universe. Since the non linear effects due to gravitational clustering occur at significantly 
smaller length scales, it is possible to use Newtonian gravity to describe these phenomena. In any region which is 
small compared to dn one can set up an unambiguous coordinate system in which the proper coordinate of a particle 
r{t) = a{t)x{t) satisfies the Newtonian equation = — Vr3> where $ is the gravitational potential. Expanding r and 
writing $ = $frw + 4> where $frw is due to the smooth (mean) density of matter and 4> is due to the perturbation 
in the density, we get 

iix + 2ax + ax = — Vr$FRW ~ ^r4> = ~Vr$FRW ~ Vx<^ (1) 

The first terms on both sides of the equation (ax and — Vr'&FRw) should match since they refer to the global expansion 
of the background FRW universe. Equating them individually gives the results 

.. ^0- ■ 1„, , la, 2ttG 9 

x + 2-x = ^V:r0 ; $FRW = -7r-»' = ^Pb?- (2) 

a a 2 a 6 

where is the gravitational potential generated by the perturbed mass density S {t, x) = [p(t^ x) — Pbit)]/ Pb(t) through 

\7l(j) = AnGa^iSp) = AirGphO^d (3) 

(The ^FRW is due to a uniform density sphere, as to be expected.) Hence, the equations for gravitational clustering 
in an expanding universe, in the Newtonian limit, can be summarized by 

Xi -I- —±, = -^Vx0; = iirGa^pbS (4) 

a a^ 

where pb{t) is the smooth background density of matter and Xi(t) is the trajectory of the i— th particle. We stress 
that, in the non-relativistic limit, the perturbed potential (f) satisfies the usual Poisson equation with the perturbed 
density contrast as the source. 

Usually one is interested in the evolution of the density contrast S{t,x.) rather than in the trajectories. Since the 
density contrast can be expressed in terms of the trajectories of the particles, it should be possible to write down 
a differential equation for (5(t, x) based on the equations for the trajectories Xi(t) derived above. It is, however, 
somewhat easier to write down an equation for 6k{t) which is the spatial Fourier transform of 5{t, x). To do this, we 
begin with the fact that the density p{x,t) due to a set of point particles, each of mass to, is given by 

p(x,t) = -^^fo[x-x,(i)] (5) 

where Xj(i) is the trajectory of the ith particle. To verify the a^'^ normalization, we can calculate the average of 
p(x, t) over a large volume V. We get 

P^it) = J — p(x,f) = ^^[y)-^-^ (6) 
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where N is the total number of particles inside the volume V and M = Nm is the mass contributed by them. Clearly 
oc a~^, as it should. The density contrast (5(x, t) is related to p(x, t) by 



l + <5(x,t)^^^ = ^^^Z3[x-x,(<)]= / rfq<5z.[x-XT(t,q)]. 
Pb N J 



(7) 



In arriving at the last equality we have taken the continuum limit by two steps: (i) We have replaced Xi(t) by XT(i, q) 
where q stands for a set of parameters (like the initial position, velocity etc.) of a particle; for simplicity, we shall take 
this to be initial position. The subscript 'T' is just to remind ourselves that XT(t, q) is the trajectory of the particle, 
(ii) We have also replaced {V/N) by d^q since both represent volume per particle. Fourier transforming both sides 
we get 



(5k(<) = J d^xc-'^ ''S{x,t) = y"d3qexp[-ik.XT(i,q)] - (27r)3,5^(k) 



(8) 



Differentiatingthis expression, and using Eq. Q for the trajectories give, after straightforward algebra, the equation 
(see 0, [nj.i): 



a a'' J ^ ' 

which can be further manipulated to give 



<5k + 2-5k = 47rGpf,(5k + - Sk 
a 



with 



Ak = 47rGpf, 



(2^ 



<^k'^k-k' 



k.k' 



(9) 



(10) 



(11) 



Bi^ — J d'^q (k.xy)^ exp [—ik.xrit, q)] 



(12) 



This equation is exact but involves XT(t,q) andxT(t,q) on the right hand side; hence it cannot be considered as 
closed. 

The structure of Eq. I|1U|) can be simplified if we use the perturbed gravitational potential (in Fourier space) (/)k 
related to 5k by 



4 = 



3H^ 



and write the integrand for Ak in the symmetrised form as 



(5k'<5k-k' 



k.k' 



= 2'^k'5k-k' 



1 fSi 



k.k' k.(k-k') 



/fc'2 |k-k'|2 



2\k'^ \ |k-k'|2 



1 / 2a 



2 V 3ffn2 



fk'fk-k' 



(k - k')2k.k' + fc'2 (fc2 _ k.k') 

fc^(k.k' + fc'2) -2(k.k')2l 



In terms of (p^ (with k' = (k/2) + p), equation H10|l becomes, 



6k + 4-(/)k 
a 



1 

'2^ 



(27r) 



gV^ik+pV'ik-p 



k^ 

k 



f3Hl\ f d^q fk.±^ ^ 



(13) 



(14) 



(15) 
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where x = xt(<, q). 

Of course, this equation is not 'closed' either. It contains the velocities of the particles iix and their positions 
explicitly in the second term on the right and one cannot — in general — express them in simple form in terms of 
0k- As a result, it might seem that we are in no better position than when we started. I will now motivate a strategy 
to tame this term in order to close this equation. This strategy depends on two features: First, extremely nonlinear 
structures do not contribute to the difference (Ak — -Bk) though, of course, they contribute individually to both 
and -Bk- Second, we can use Zeldovich approximation to evaluate this term, once the above fact is realised. I will now 
elaborate on these two features. 



A. 'Renormalizability' of gravity 

Gravitational clustering in an expanding universe brings out an interesting feature about gravity which can be 
described along the following lines. Let us consider a large number of particles which are interacting via gravity in 
an expanding background and forming bound gravitating systems. At some time t, let us assume that a fraction / 
of the particles are in virialized, self-gravitating clusters (of typical size R, say) which are reasonably immune to the 
effect of expansion. Imagine that we replace each cluster by a single particle at its center of mass with the mass 
equal to the total mass of the cluster. (The total number of particles have now been reduced but, if the original 
number was sufhciently large, we may assume that the resulting number of particles is again large enough to carry on 
further evolution with a valid statistical description.) We now evolve the resulting system to a time t' and compare 
the result with what would have been obtained if we had evolved the original system directly to t' . Obviously, the 
characteristics of the system at small scales (corresponding to the typical size R of the clusters at time t) will be quite 
different. However, at large scales (kR -C 1), the characteristics will be the same both the systems. In other words, 
the effect of a bunch of particles, in a virialized cluster, on the rest of the system is described, to the lowest order, by 
just the monopole moment of the cluster - which is taken into account by replacing the cluster by a single particle 
at the center of mass having appropriate mass. In this sense, gravitational interactions are "renormalizable" - where 
the term is used in the specific sense defined above. 

The result has been explicitly verified in simulations but one must emphasize that the whole idea of numerical 
simulations of such systems tacitly assumes the validity of this result. If the detailed non linear behaviour at small 
scales, say within galaxies, influences very large scale behaviour of the universe (say, at super cluster scales), then it 
will be impossible to study the large scale structure in the universe with simulations of finite resolution. 

One may wonder how this feature (renormalizability of gravity) is taken care of in Eq. (|10|l . Inside a galaxy cluster, 
for example, the velocities can be quite high and one might think that this could influence the evolution of (5k at 
all scales. This does not happen and, to the lowest order, the contribution from virialized bound clusters cancel in 
— ^k- We shall now provide a proof of this result (also see 0|). 

We begin by writing the right hand side TZ of the Eq. @ concentrating on the particles in a given cluster. 



TZ = J d3qe-* Q,r |^+zQ'^(k-Q)| 



(16) 



where we have used the notation Q = xy for the trajectories of the particles and the subscripts a, b, .... = 1, 2, 3 denote 
the components of the vector. For a set of particles which form a bound virialized cluster, we have from Eq. gl) the 
equation of motion 

Q' + 2«Q'.-4|^ (17) 

We multiply this equation by Q-', sum over the particles in the particular cluster and symmetrize on i and j, to obtain 
the equation 

We use the summation symbol, rather than integration over q merely to emphasize the fact that the sum is over 
particles of a given cluster. Let us now consider the first term in the right hand side of Eq. IjlGI) with the origin 
of the coordinate system shifted to the center of mass of the cluster. Expanding the exponential as e"''' *^ w 
(1 — ik - Q) + 0{k^R?) where R is the size of the cluster, we find that in the first term, proportional to V0, the 
sum of the forces acting on all the particles in the cluster (due to self gravity) vanishes. The second term gives, on 
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symmetrization, 
Using Eq. (^HJ wc find that 

^a-(k.V0)e--Q = + 5:(k . Q)^ + + 2^1) • Q)^ (20) 

The second term is of order 0{k^R^) and can be ignored, giving 

ia-\-k-V^)e-'^-^ « + ^(k • Qf + 0{k^R^) (21) 

Consider next the second term in the right hand side of with the same expansion for the exponential. We get 

J2{ik'')e-'^-^ [iQ"fc''Q^j = - E k''k''Q''Q''{l - ik ■ Q) + 0{k^R^) 

= - 2(k • Q)2 + ^(k ■ Q)2fc'^Q'^ + OiPR^) (22) 

The second term is effectively zero for any cluster of particles for which Q ^ — Q is a symmetry. Hence the two 
terms on the right hand side of Eq. pt)|) cancel each other for all particles in the same virialized cluster; that is, to 
the order 0{k^R^), the term (Ak — -Bk) receives contribution only from particles which are not bound to any of the 
clusters. If the typical size of the clusters formed at time t is R, then for wave-numbers with k^R^ <^ 1, we can ignore 
the contribution from the clusters. Hence, in the limit of fc ^ we can ignore (Ak — -Bk) term and treat equation 
Eq. as linear in (5k; large spatial scales in the universe can be described by linear perturbation theory even when 
small spatial scales are highly non linear. 

There is, however, an important caveat to this claim. In the right hand side of Eq. (|10|l one is comparing the first 
term (which is linear in 5k) with the contribution (Ak — Bu). If, at the relevant wavenumber, the first term 47rGpf,^k 
is negligibly small, then the only contribution will come from (A^ — B^) and, of course, we cannot ignore it in this 
case. The above discussion shows that this contribution will scale as k^R^ and will lead to a development of oc fc^ 
if originally (in linear theory) 6k oc fc" with n > 2 as A: — > 0. We shall say more about this later on. 

Clearly, we can ignore the contribution from particles in virialized clusters in estimating (^k ^ -Bk). We will next 
consider how we can estimate the effect of remaining particles. A useful insight can be obtained by examining the 
nature of particle trajectories which lead to the growth of the density contrast (5k oc a in the linear limit. To determine 
the particle trajectories corresponding to the linear limit, let us start by writing the trajectories in the form 

XT(a,q) = q + L(a,q) (23) 

where q is the Lagrangian coordinate (indicating the original position of the particle) and L(a, q) is the displacement. 
The corresponding Fourier transform of the density contrast is given by the general expression 

5k(«) = J d3qe-* 'i-^'' L(-.q) - (2^)3(5Di,,,[k] (24) 

In the linear regime, we expect the particles to have moved very little and hence we can expand the integrand in the 
above equation in a Taylor series in (k • L). This gives, to the lowest order, 

dkia) = -J d^qe-*i(ik-L(a,q)) = d^qg-^ q ■ L) (25) 

showing that (5k(a) is Fourier transform of — Vq.L(a, q). This allows us to identify V • L(a, q) with the original density 
contrast in real space —Sq{a). Using the Poisson equation we can write Sq{a) as a divergence; that is 

V • L(a, q) = -Sq{a) = - ^^o^'aV • ( V(/)) (26) 

which, in turn, shows that a consistent set of displacements that will lead to (5(a) oc a is given by 

L(a, q) = -(V^)a = au(q); V = (2/3)i?o~'<^ (27) 
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The trajectories in this hmit are, therefore, hnear in a: 

XT (a, q) = q + au(q) (28) 

A useful approximation to describe the quasi hnear stages of clustering is obtained by using the trajectory in Eq.J^HJ 
as an ansatz valid even at quasi linear epochs. In this approximation, (called Zeldovich approximation), the proper 
Eulerian position r of a particle is related to its Lagrangian position q by 

v{t) = a{t)x{t) = a{t) [q + a(t)u(q)] (29) 

where x(i) is the comoving Eulerian coordinate. If the initial, unperturbed, density is p (which is independent of q), 
then the conservation of mass implies that the perturbed density will be p{r,t)d'^r = pd^q. As is well known, this 
suggests that sheet like structures, or 'pancakes', will be the first nonlinear structures to form when gravitational 
instability amplifies density perturbations. 



B. Closure ansatz for the dynamical equation 



We now combine the two results obtained above in order to suggest a closure condition for our dynamical equation. 
We begin by noting that at any given moment of time we can divide the particles in the system into two sets - 
those which are already a part of virialized cluster and those which are not. Of these, we know that the first set of 
particles do not contribute significantly to {A^ — i^k) so we will not incur any serious error in ignoring these particles 
in computing (Ak — i?k). For the description of particles in the second set, the Zeldovich approximation should be 
fairly good. In fact, we can do slightly better than the standard Zeldovich approximation. We note that in Eq. H29|l 
the velocities were taken to be proportional to the gradient of the initial gravitational potential. We can improve on 
this ansatz by taking the velocities to be given by the gradient of the instantaneous gravitational potential which has 
the effect of incorporating the influence of particles in bound clusters on the rest of the particles to certain extent. 
Given this ansatz, it is straightforward to obtain a closed integro-differential equation for the gravitational potential 
along the following lines. The trajectories in Zeldovich approximation, given by Eq. (EHJ leads to: 

XT(q,a) = q + aVV'; xt = V?/'; V' = g^"^ (30) 

To the same order of accuracy, i?k in Eq. H12|l becomes: 

J d\ (k • Xt)^ e-* ('i+^) ^ J d^q{k ■ ±T)'^e-'^ '^ (31) 

Substituting these expressions in Eq. ((T^ we find that the gravitational potential is described by the closed integral 
equation: 

,< _ _J_ f d^p 
a"^" ~ 3a2 7 (2^)3 



/.k + 4-,^k = I 7:^01 k+p^ik-p^lk,?) 



g(k,p) = le + ^-p^-B^^-"" 
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(32) 



This equation provides a powerful method for analyzing non linear clustering since estimating (Ak — -Bk) by Zeldovich 
approximation has a very large domain of applicability. 

At this stage, the validity of the above equation rests on the conjectured validity of the ansatz used in approximating 
the trajectories. Though I have given a brief argument to motivate this ansatz, its ultimate validity can be tested only 
by numerical simulations as well as by working out the consequences of the equations. The rest of the paper works out 
several analytical results which can be obtained from the above equation. We have also performed detailed numerical 
simulations to check the above ansatz and some of these results will be presented elsewhere Q. For the propose of 
this paper, I merely state that the numerical simulations show that this ansatz has a wide domain of applicability. 

In the next two sections, I will use this equation to study the following important questions: Suppose power was 
injected into a self gravitating system at a fixed length scale [that is, the initial power spectrum is Pin(k) oc (5D(|k|— fco)]. 
How does the non linear evolution transfer the power to other wave numbers? The corresponding question in fluid 
turbulence is very well studied and we know that for a wide range of scales the final spectrum is the standard 
Kolmogorov spectrum. In the next two sections, we shall investigate some aspects of this problem using the above 
formalism. 
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III. INVERSE CASCADE IN NON LINEAR GRAVITATIONAL CLUSTERING: THE fc* TAIL 



There is an interesting and curious result which is characteristic of gravitational clustering that can be obtained 
directly from our Eq. H15|l . Consider an initial power spectrum which has very little power at large scales; more 
precisely, we shall assume that P{k) dies faster than fc^ for small k. If these large scales are described by linear theory 
— as one would have normally expected — then the power at these scales can only grow as and it will always be 
sub dominant to fc^. It turns out that this conclusion is incorrect. As the system evolves, small scale nonlinearities 
will develop in the system and — if the large scales have too little power intrinsically (i.e. if n is large) — then the 
long wavelength power will soon be dominated by the "tail" of the short wavelength power arising from the nonlinear 
clustering. This occurs because, in Eq. (|10|l . the nonlinear term (Ak - Bk) = 0{k^R^) can dominate over AnGpbSk 
at long wavelengths (as k ^ 0) and lead to the development of a k'^ power spectrum at large scales. This is a purely 
non linear effect which we shall now describe. 

A formal way of obtaining the /c^ tail is to solve Eq. (|32|l for long wavelengths; i.e. near k = 0. Writing (py^ = 

"^^k ' + + .... where 0^^^ = (f)^"^ is the time independent gravitational potential in the linear theory and (/)^^^ is the 
next order correction, we get from Eq. H32() . the equation 

The solution to this equation is the sum of a solution to the homogeneous part [which decays as </> oc a^^ oc t^^^^ 
giving (/) oc and a particular solution which grows as a. Ignoring the decaying mode at late times and taking 

'^k ~ one can determine Ck from the above equation. Plugging it back, we find the lowest order correction to 
be, 

*^"^-(2ij)/|^Wt.-.«..P) (34, 

Near k ~ 0, we have 

,(2) _ 2a f d^p ^ 2 

OO 



(35) 

which is independent of k to the lowest order. Correspondingly the power spectrum for density Ps{k) oc a?'k'^P^{k) oc 
a^fc^ in this order. 

The generation of long wavelength k'^ tail is easily seen in simulations if one starts with a power spectrum that 
is sharply peaked in |k|. Figure shows the results of such a simulation (adapted from @) in which the y-axis is 
[A(fc)/a(i)] where A^(fc) = k^P/2Ti^ is the power per logarithmic band in k. In linear theory A oc a and this quantity 
should not change. The curves labelled by a = 0.12 to a = 20.0 show the effects of nonlinear evolution, especially the 
development of fc'* tail. 



5(k.p)^ 



IV. CASCADING IN NON LINEAR GRAVITATIONAL CLUSTERING: GENERATION OF SMALL 

SCALE POWER 

We next turn to the transfer of power to smaller spatial scales due to nonlinear mode coupling. Figure ^ also shows 
that, as the clustering proceeds, power is generated at spatial scales smaller than the scale k^^ at which the power is 
injected. (One can also see a characteristic peak at smaller scales.) These features can also be easily understood from 
Eq. H34I) . Let the initial gravitational potential and the density contrast (in the linear theory) be sharply peaked at 
the wave number fco, say, with: 

0k = I'^^dM - ^o]; kl5^ = -^{txa)ko5Dm - ^o] (36) 

where /x is dimensionless constant indicating the strength of the potential and the other factors ensure the correct 
dimensions. Equation H34|l shows that, the right hand side is nonzero only when the magnitudes of both the vectors 
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FIG. 1: The transfer of power to long wavelengths forming a tail is illustrated using simulation results. Power is injected 
in the form of a narrow peak at L = 8. Note that the j/— axis is (A/a) so that there will be no change of shape of the power 
spectrum under linear evolution with A oc a. As time goes on a k* tail is generated purely due to nonlinear coupling between 
the modes. (Figure adapted from ref.0.) 



[(l/2)k + p] and [(l/2)k — p] are fco- This requires k-p = and (fc/2)^+p^ = fcg. (Incidentally, this constraint has a 
simple geometric interpretation: Given any k, with k < 2fco one constructs a vector k/2 inside a sphere of radius fco 
and a vector p perpendicular to k/2 reaching up to the shell at radius fcg where the initial power resides. Obviously, 
this construction is possible only for k < 2ko.) Performing the integration in Eq. H34() we find that 



(We have again ignored the decaying mode which arises as a solution to the homogeneous part.) The corresponding 
power spectrum for the density field P(fc) = |<5/cp oc a^fc*|0fcp will evolve as 

p(2)(fc)a(HVfi-^9')'(i + ig')'; g = ^ (38) 



The power at large spatial scales (fc 0) varies as fc** as discussed before. The power has also been generated 
at smaller scales in the range fco < fc < 2fco with P*^^)(fc) being a maximum at k„i ~ 1.54fco corresponding to the 
length scale k^^ w 0.65fc(^^. Figure |21 shows the power spectrum for density field (divided by to eliminate linear 
growth) computed analytically for a narrow Gaussian initial power spectrum centered at fcp = 1. The curves are for 
(/ia/567r^)^ — 10~^, 10~^, 10~^ and 1. The similarity between figures ^ and |21 is striking and allows us to understand 
the simulation results. The key difference is that, in the simulations, newly generated power will further produce power 
at 4fco, 8ko, ... and each of these will give rise to a fc^ tail to the right. The resultant power will, of course, be more 
complicated than predicted by our analytic model. The generation of power near this maximum at fc~^ = 0.65kQ^ is 
visible as a second peak in figure El and around 27r/fco « 4 in figure^ 

If we had taken the initial power spectrum to be Dirac delta function in the wave vector k (rather than on the 
magnitude of the wave vector, as we have done) the right hand side of Ea. (|34|l will contribute only when (^kztp) = ko. 
This requires p = and k = 2ko showing that the power is generated exactly at the second harmonic of the wave 
number. Spreading the initial power on a shell of radius fcoj spreads the power over different vectors leading to the 
result obtained above. 
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FIG. 2: Analytic model for transfer of power in gravitational clustering. The initial power was injected at the wave number 
ko with a Gaussian window of width Afc/fco = 0.1. First order calculation shows that the power is transfered to larger spatial 
scales with a fc* tail and to the shorter spatial scales, all the way down to (l/2)fc(^^. The plot gives the total power spectrum 
divided by (with y— axis normalized arbitrarily) at different times with changing by factor 10 between any two curves. 

Equation H3t)|l shows that fcgfJk will reach non linearity for w (3/2). The situation is different as regards the 
gravitational potential due to the large numerical factor SGtt^; the gravitational potential fluctuations are comparable 
to the original fluctuations only when /ia « 56tt^. 

V. CONCLUSIONS 

The purpose of this paper was to draw attention to a possible approach to study nonlinear gravitational clustering. 
I showed how one can obtain a closed integro-differential equation for the evolution of the gravitational potential in the 
Fourier space. This equation shows that the nonlinear evolution of the potential at the wave number k is essentially 
dictated by a two-mode coupling between modes at wave numbers (l/2)k + p and (l/2)k — p integrated over all other 
modes p with a quadratic kernel. Among other things, this equation is useful in studying how the power injected at 
a given scale flows to other scales in gravitational clustering. It is seen that the nonlinear evolution leads to a k'^ tail 
at larger spatial scales and pumps power into smaller scales through repeated generation of higher harmonics. 

The derivation of the dynamical equation was based on an ansatz which needs to be verified by simulations. We 
shall report on this aspect as well as on more detailed features of power transfer in a separate publication 9] . There are 
several other obvious directions in which this formalism can be developed further. For example, one can incorporate 
the effect of bound clusters by adding the (pi,: oc k^^ term due to the monopole moments. One can also improve the 
accuracy by separating the contribution from different scales in the Fourier space into nonlinear (k > fc„/(t)) and 
linear (fc < kni{t)) scales (where kni{t) is the scale that is going nonlinear at time t) and dealing with them separately. 
These issues are under investigation. 
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APPENDIX A: SPHERICAL APPROXIMATION 



The purpose of this brief appendix is to show how the formahsm developed in the text connects up with standard 
spherical approximation used very often in cosmology. In our language, this approximation consists of assuming that 
the trajectories are homogeneous; i.e. x(t, q) = /(i)q where f{t) is to be determined. In this case, the density contrast 
is 

<5k(i) = J d^qe-^/Wk-q - {2nf5D{k) = i2nf Somr' " 1] = i2Trf 60(^)5 (t) (Al) 

where we have defined S{t) = [f^^{t) — l] as the amplitude of the density contrast for the k = mode. It is now 
straightforward to compute A and B in Ea. (|l()(l . We have 



and 



A^A7TGpbS^{t)[i27Tf5Dm (A2) 



B = jd\{k-q,fpe-^f'^^^<i^) = -p-^[{2^f5D{m 

^ [{2i.f5nm (A3) 



so that the Eq. (|10|l becomes 



3(1 + ^)' 



^ + 2-(j = 47rGpfe(l + (5)<5+|-^ (A4) 
a 3 (1 + dj 



It is easy to show that this equation is identical to that of standard spherical approximation governed by i? = ~GM/R^ 
if we take (1 + (5) oc (a/R)^ . Thus our formalism allows one to reproduce all the known, standard approximations for 
nonlinear epochs. 
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